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We extend the variational cluster approach to deal with strongly correlated lattice bosons in the 
superfluid phase. To this end, we reformulate the approach within a pseudoparticle formalism, 
whereby cluster excitations are described by particlelike excitations. The approximation amounts 
to solving a multicomponent noninteracting bosonic system by means of a multimode Bogoliubov 
approximation. A source-and-drain term is introduced in order to break U(l) symmetry at the 
cluster level. We provide an expression for the grand potential, the single-particle normal and 
anomalous Green's functions, the condensate density, and other static quantities. As a first nontrivial 
application of the method we choose the two-dimensional Bose-Hubbard model and evaluate results 
in both the Mott and the superfluid phases. Our results show an excellent agreement with quantum 
Monte Carlo calculations. 
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I. INTRODUCTION 

Cluster approaches have been proven to be very use- 
ful for the numerical investigation of strongly correlated 
many-body systems. These approaches consist in em- 
bedding finite size clusters, for which a numerical exact 
solution is available, within a lattice of infinite size. The 
embedding is done by introducing additional fields to the 
cluster Hamiltonian, in order to take into account the 
coupling to the rest of the lattice in some appropriate dy- 
namical mean-field way. We will term these fields Weiss 
fields, since they play an analogous role as in Weiss mean- 
field theory of ferromagnetism (c/. Ref.[l|). Different clus- 
ter embedding techniques, such as cluster perturbation 
theory^ (CPT), variational cluster approach^— (VCA) 
, cellular dynamical mean field theory^ (C-DMFT), and 
dynamical cluster approximation^ differ by the nature 
of the Weiss fields and of the mean-field treatment which 
fixes their optimal value. In the present paper, we con- 
sider VCA, which has been applied to a large variety of 
fermionic&£ and bosonio^r— systems. VCA can be un- 
derstood in a more general framework called self-energy 
functional approach 5 *^ (SFA), in which the grand po- 
tential of the physical system is expressed as the station- 
ary point of a particular functional of the self energy. 
Here, we will adopt an alternative approach to VCA in 
which single-particle excitations are expressed in terms 
of "pseudoparticlcs," which arc similar to Hubbard op- 
erators^ and external fields are "added" to the clus- 
ter Hamiltonian and "subtracted perturbatively."— We 
discuss in Sec. |H] the advantages of this alternative ap- 
proach. 

Strongly correlated lattice bosons are currently in the 
focus of research due to seminal experiments on ultracold 
gases of atomsJ^— In these experiments quantum me- 
chanical interference effects can be observed on a macro- 
scopic scale. In particular, ultracold gases of atoms on 
a lattice undergo a quantum phase transition from the 



Mott phase, in which particles are localized on individ- 
ual lattice sites, to the delocalized superfluid phase, in 
which t/(l) symmetry is broken and a finite fraction of 
the particles forms a Bose-Einstein condensate (BEC). 

Up to now, bosonic VCA has been formulated for the 
normal phase only. The principal aim of this paper is 
to extend this formalism to the symmetry-broken, super- 
fluid phase. The theoretical framework, developed in the 
next two sections, is applicable to a large class of lat- 
tice boson systems in the Mott insulating as well as in 
the superfluid phase. In particular, besides the widely 
studied Bose-Hubbard model ) 15 ' 18 the method can be 
straightforwardly extended to include disordered systems 
or multiple components containing, for example, fcrmion- 
boson mixtures. The extended VCA theory can be ap- 
plied even to the U(l) broken, superfluid phase of light- 
matter systems, where photons are confined in coupled, 
nonlinear quantum-electrodynamics cavitiesi 19 : 20 In or- 
der to achieve the extension to the U(l) broken, su- 
perfluid phase, it proves convenient to reformulate VCA 
in terms of a pseudoparticle approach, whereby single- 
particle excitations within a cluster arc approximately 
mapped onto particlelike excitations. We show that this 
approach, first applied to normal bosons, quite naturally 
suggests the extension to the superfluid case. In a follow- 
ing publication^ we show that the results obtained from 
the pseudoparticle formalism in the superfluid phase can 
be equivalently obtained within an appropriate extension 
of the SFA taking into account condensed bosons. One 
of the aims of the present paper is to illustrate the ad- 
vantages of the pseudoparticle formalism, which can be 
used to extend VCA to a large variety of problems with 
strongly correlated lattice systems. 

The pseudoparticle formalism is in some aspects re- 
lated to the standard basis matrix operator method 
developed by Haley and Erdos in Ref. [22| and to the 
Hubbard-operator approach, sec for instance Ref. [l3|. 
The idea is to introduce pseudoparticle operators and 
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&L which connect the ground state |?/>o) with single- 
particle excited states |^>„) of a Hamiltonian describing 
disconnected clusters in the lattice. In the VCA language 
the cluster Hamiltonian is termed reference system H'. 
Of course, the 6 M and 6J, do not have the properties of 
ordinary single-particle creation and annihilation opera- 
tors. The crucial point is that by treating them as such, 
one recovers the very same results as obtained from CPT 
and from VCA, as has been shown for fermionic systems 
in Ref. [HI (see appendix therein). In Sec. [H]we prove the 
same result for the bosonic (normal) case, which is some- 
what more subtle, as it requires a multimode Bogoliubov 
transformation. In this picture, excited states are 
treated as pseudoparticle excitations with the properties 

|W) = bl |Vo> = <W ■ 

Within the VCA approximation, pscudoparticles are re- 
garded as nonintcracting particles. We stress that, while 
this may seem a rather crude approximation, it is equiv- 
alent to CPT and VCA. Furthermore, with appropriate 
extensions it becomes equivalent to C-DMFT. 

It is straightforward to show (see Sec. [TTJ) that within 
this approach the original bosonic operators a i and aj 
can be expressed as linear combinations of the pseu- 
doparticle operators b^. This makes it possible to write 
the coupling of the cluster to the rest of the lattice, 
which in VCA consists of intercluster hopping terms, as 
a quadratic form in the b^. In combination with the fact 
that the cluster Hamiltonian is by construction quadratic 
in these operators as well, one finally obtains a Hamilto- 
nian which is completely quadratic in the pseudoparticle 
operators, and can, thus, be solved exactly. 

Our paper is organized as follows: In Sec. |H] we first 
show that the standard VCA results for the Green's func- 
tions and for the grand potential fi are recovered within 
the pseudoparticle formalism applied to normal bosonic 
systems. In order to be able to treat the supcrfluid 
phase we then extend the theory in Sec. Mil by introduc- 
ing Weiss fields in the form of "source-and-drain" terms 
which explicitly break the U(l) symmetry of the refer- 
ence Hamiltonian H'. The main result of Sec. IIIII is the 
expression for the grand potential N c fl of the physical 
system ( N c is the number of clusters) . Within our exten- 
sion of VCA to the supcrfluid phase we obtain^ 

Q = n' — Trln(-G) + — Trln(-G') - -tih 

2N C K ' 2N C V ! 2 

+ 1 -{A^)G^ ) {A)- 1 -{A^G'^{A)' . (1) 

The first three terms on the right hand side are essen- 
tially identical to those which arc also present in standard 
VCA expressions. Particularly, 57' is the grand potential 
(per cluster) of the reference system, and G' and G are 
the connected Green's functions of the reference and of 
the physical system, respectively. However, they are ex- 
pressed in the Nambu representation, which explains the 
additional factor 1/2 and the fourth term in comparison 



with previous results^— The suffix (0) used in the sec- 
ond line of Eq. ([T]) means that the corresponding Green's 
functions are calculated for q = and w = 0, where q 
is the superlattice vector associated to the cluster tiling, 
and ui is the Matsubara frequency. As usual within VCA 
theory, the two Green's functions share the same self- 
energy. The expectation values (A)' and (A) are the cor- 
responding condensate densities, again in Nambu (vec- 
tor) notation. The latter are connected by the relation 

G^(A) = F + G{-i(Ay , (2) 

where the vector F describes the strength of the source- 
and-drain term which is introduced in the reference sys- 
tem in order to explicitly break U(l) symmetry. The 
value of F [see Eq. ([23)) ] has to be determined from the 
variational principle. Details for the notation are pro- 
vided in Sec. [Hi and IIIII In addition to the formula 
for the grand potential VL, we evaluate expressions for 
other quantities, which are useful for describing the su- 
pcrfluid phase. In particular, wc derive expressions for 
the normal and anomalous Green's functions, the parti- 
cle density, and the condensate density. In Sec. IIVI this 
extended VCA theory is applied to the two-dimensional 
Bose-Hubbard (BH) model in the superfluid phase. Fi- 
nally, we summarize and conclude our findings in SecfVl 



II. PSEUDOPARTICLE APPROACH 

In this section we reformulate CPT/ VCA within the 
pseudoparticle approach for bosonic systems. In princi- 
ple, one may argue that the formulation of CPT/ VCA 
using pscudoparticles is complicated and in the case of 
the normal phase (i.e. Mott phase) CPT/ VCA can be 
obtained from simpler approaches, as, for example, from 
Dyson's equation (see, e.g. Ref. 0) and the SFA££ The 
reason why we present this alternative formulation here 
is that this approach, while not as rigorous as SEA, pro- 
vides useful hints on how to deal with more complicated 
situations, like the supcrfluid phase discussed in this work 
(see Sec. IIII|). In addition, it gives insight on other prop- 
erties. For example, in the case of normal bosons the 
pseudoparticle approach is useful in order to understand 
the occurrence of noncausality of the Green's function in 
cases, where the chosen reference system is not suitable 
to describe the phase of the physical system, as we point 
out below. Thus the aim of this section is to derive the 
principal theoretical framework of the pseudoparticle ap- 
proach, for the normal phase, reproducing the known re- 
sult for the grand potential fl, which has to be optimized. 
The extension to the superfluid phase is the subject of 
the next section. 

The physical system of interacting particles is de- 
scribed by a grand-canonical Hamiltonian H, which is 
related to the canonical Hamiltonian in the usual way by 
the additional single-particle term —fiN. The Hamiltoni- 
ans, which can be treated by the extended VCA theory, 
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generally have the form H = H t + Hjj, where H t con- 
sists of arbitrary one-particle terms and Hjj of local two- 
particle terms. The physical system is defined on a large 
or even infinite lattice with periodic boundary conditions. 
The underlying lattice is now tiled into N c clusters each 
one containing L orbitals (sites). We split the Hamilto- 
nian into a cluster part H c i, which only describes pro- 
cesses within the various clusters, and the residual part 
T, containing the intercluster processes, which consist of 
single-particle terms only, so that 

H = H cl +f. (3) 

CPT amounts to first solving for the Hamiltonian H c i 
and then carrying out a perturbation expansion in the 
intercluster Hamiltonian T . Of course, within CPT one 
is free to add an arbitrary single-particle Hamiltonian 
—A to the cluster Hamiltonian H c i provided it is then 
subtracted from T so that H remains unchanged. This 
defines a new cluster Hamiltonian H' 

H' = Hd - A . (4) 

The physical Hamiltonian H , given in Eq. ([3]), can now 
be expressed in terms of the new cluster Hamiltonian 

H = H' + A + f = H' + f, (5) 

leading to a new "perturbation" f = A + f . The CPT 
expansion is now carried out in this new "perturbation" . 
While ideal exact results should not depend on A (this 
occurs, for example, in the noninteracting case), in prac- 
tice results do depend on A due to the approximate na- 
ture of the expansion. The idea is to fix the parame- 
ters A by an optimization prescription, which amounts 
to finding the stationary point of the grand potential 
obtained from the perturbative expansion. The optimiza- 
tion prescription is put on a rigorous framework within 
the SEA.— It is straightforward to show that this pro- 
cedure is equivalent to the standard VCA prescription, 
whereby H 1 is the corresponding reference system.— 

In the following, we consider N c identical disconnected 
clusters, and denote the sites (orbitals) within a cluster 
by i. The position of each cluster on the large, physical 
lattice is specified by a lattice vector R. Accordingly, we 
denote by o^r the annihilation operator for a boson on 
site i of cluster R, and similarly for creation operators 
a \ r- I R 01 'der to keep a compact notation we combine the 
annihilation operators of a given cluster R into a column 
vector of operators 

aR = (ffli,Rj <J2,R, • ■ • a L,n) T , 

and correspondingly, the creation operators are row vec- 
tors a R = (aa)t Using these expressions we rewrite the 
intercluster Hamiltonian as 

f= ]Ta R i(R-R)a R , , (6) 

RR' 



where i(R— R') is a matrix describing the hopping terms 
from cluster R' to cluster R, with the property t(R — 
R') = i(R' — R)' . Here we have assumed translation 
invariance by a cluster translation vector. Similarly, we 
can express A in terms of an intracluster hopping matrix 
h 

R 

such that T, defined in Eq. ([5]), can be written as Eq. (|6]) 
with the replacement 

t(R - R) -> i(R - R) = t(R - R) + 5 R . R ,h . 

As explained above, the reference system consists of a 
sum of Hamiltonians acting on independent clusters R 

R 

Again considering translation invariance, all H'(TV) are 
identical. Thus it suffices to determine numerically the 
ground state |V>0jR)j a s well as single particle or single- 
hole excited states | ?/> M , R) of a single cluster Hamiltonian 

ir(R), with corresponding eigenenergies E' Q and E' re- 
spectively. The key idea of the approach, to be presented 
here, is to introduce pseudoparticle operators b^ R and 

b R = (&^ R )^, which are defined by their matrix ele- 
ments 

(W,R|&£, R hfo.B.} = V- (?) 

In other words, the pseudoparticle operator b^ R applied 
to the exact many-body groundstatc |?/»o,R) of a clus- 
ter creates the exact excited many-body state |i/) ([i ,R). 
In this respect, it is of course forbidden to apply a 
second pseudoparticle creation operator on the excited 
state. This leads to the supplementary hard-core con- 
straints b v R &|j R IV'Oj R) =0. To neglect this hard-core 
constraint and to restrict to single- particle and single- 
hole excitations within each cluster is the approxima- 
tion made here. We show below that this approxima- 
tion, combined with the variational procedure discussed 
above, gives the same results as VCA. In particular, we 
obtain the same expression for the grand potential f2, and 
for the Green's function. It should be mentioned, how- 
ever, that within the pseudoparticle approach there is no 
known rigorous variational principle for f2. One can sim- 
ply heuristically state that the "best" solution is the one 
that "minimizes" the energy, although, as we know from 
VCA, the variational solution is not always a minimum. 
Also for parameters, such as the chemical potential, for 
which f2 turns out to be a maximum, one can argue that 
the stationary condition is a kind of "constraint" fixing 
the consistency of thermodynamic quantities,— and the 
corresponding parameter is a kind of "Lagrange multi- 
plier." Nevertheless, it is not the goal of the present pa- 
per to discuss this issue. Here, we want simply use this 
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"tool" in order to formulate an extension of the theory to 
address the bosonic superfluid phase (see Sec. IIIip . The 
knowledge of the correction to the order parameter and 
of the grand-potential f2 can then guide and facilitate 
a rigorous extension of SFA to deal with the superfluid 
phase. This is the goal of a future publication.— 

With the help of these operators, it is straightforward 
to write down a Hamiltonian which has the same ener- 
gies and eigenvectors as the reference system, restricted 
within the subspace of single-particle and single-hole ex- 
citations from the ground state 

H' = N c Sl' + Y^Y, AE '» b UK,iv (8) 

R v 

with the (positive) excitation energies AE' V = E' v — E' Q . 
Since we are interested in zero temperature T = 0, the 
grand potential of the reference system is Q' = E' . 

To proceed further, we need an expression for T, and, 
thus, of the original bosonic operators a^R,, in terms of 
the pseudoparticle operators. For simplicity, we drop the 
R index and concentrate on a given cluster. Within the 
pseudoparticle approximation the operators must coin- 
cide only within the constrained subspace. We thus ap- 
proximate each a,i by an operator O i (b^, bjj) which shares 
the same matrix elements (-00 1 ■ |Vo/i (Vol ' IVf)j and 
(Vv| • |Vo)- We express Oi by means of the ansatz 



0<(mJ.) = XXa+ 



E 

fi=n p + l 



where the first sum contains the n p indices associated 
with the single-particle excitations, and the second sum 
contains the indices for the single-hole excitations. 
The total number of excitations taken into account is 
n s = n p + nh- Here we have exploited particle-number 
conservation. Next, we use this expression to evaluate 
the following matrix elements 



(Vol 6, (V&t^oH^ =(ih\ai\fM (10a) 
(Vvl O i (&„, 6+) |Vo) = Zi, v = m a t |Vo) (10b) 
(Vo| O t (b M , b\) |V,) = Ri,» = (Vol <H |V.) , (10c) 

where the coefficients r y i are zero so far, since the ref- 
erence system conserves the particle number. We now 
introduce the compact notation 



t 



B = (b 1> ...,b np> b{, 1 ...bi) T & = {B) 



i.e., the first part of the vector acts on particle states, 
and the second part on hole states. Notice that in this 
form (B) changes the number of particles by +1 (—1). 
We also introduce the Q matrix (which is the same as in 
Ref. [ll|) as 

/ Ri,v for 1 < v < n p 
7 ji v for n„ < v < n s 



The Q matrix can be used to express the original oper- 
ators a and a' in terms of B operators [cf. Eq. ©] in a 
compact form: 



a =QB 
a 1 = fitQt 



(11a) 
(lib) 



Using the compact vector notation for B and B\ the 
reference Hamiltonian [Eq. (JSJ] can be written as 

H' = N C Q' + B^SABr - N c AE' h , (12) 

R 

where we reintroduced the R dependence. Here we in- 
troduced the diagonal matrices 

S = diag(l, ... ,1, -1, . -1 ) 

1, ... ,n p rip+1, ... .n s 

and 

A = 5diag(A^, . . . , AE' np , AE' np+1 ...,AE' n J. 

Notice that S 2 — 1, while A contains the poles of the 
Green's function for the reference system. The constant 



+ % 1 , (9) w ith the function 



E 

fi=n p + l 



g(e) = ee(-e) 



tr«?(A) , 



takes into account that some of the boson operators have 
been rearranged in order to obtain Eq. (fT2j) . The physical 
Hamiltonian introduced in Eq. ([5]) reads 

R,R' 

Using Eqs. (jTTJ) and (fT2|) yields a quadratic expression in 
the B operators: 

H = N C Q' + N c tr g(A) + ^ B^SABr 

R 

+ E B * & - R ') Q Bri ■ 

R,R' 

We can now introduce a Fourier transform in the cluster 
vectors R 



(&l,q, . . . , &n p ,q, b\ p + 1 _ q , . . . , & nai _ q ) 



leading to 



H = N C Q' + N c tr g(A) + J2H q 



(13) 



(14) 
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with 

P q EE P q SM q P q . (15) 

Here, we have introduced the matrix 

M q ee A + S Q f i q Q , 

where 

R 

is the Fourier transform of t(R— R'). The non-Hermitian 
matrix M q is identically defined as in Ref. [111. 

Being quadratic in the B operators, P q can be quite 
generally put into diagonal form by a multimode Bo- 
goliubov transformation. To achieve this, we look for 
"normal-mode" pscudoparticles described by the vector 
P with the same structure as P (in the following consid- 
erations we omit the q dependence for simplicity) 

P* = (p(i,...p<) P = (pt)t, 

where s[ = ±1 so that pf 1 = p\ is a creation and pj 1 = pi 
is an annihilation operator. The new P operator shall be 
connected with B via 

B = VP , 

where V is a nonsingular but in general nonunitary ma- 
trix. From a physical viewpoint the nonsingularity of 
V corresponds to a pseudoparticle conservation, mean- 
ing that there are as many pseudoparticles B as normal- 
mode pseudoparticles P. The transformation V must 
satisfy two conditions. First it must be chosen such that 
P has appropriate bosonic commutation relations, i. e., 

[ppt] = S'EEdiag( S i,...,0- 

This gives 

s"=[p,pt] = v- 1 [b,b^](v- 1 )^ = v^siv-y , 

which in turn yields 

V S'V* S = I (16a) 
5" F f S = V~ x (16b) 
V ] SV = S' . (16c) 

The second requirement on V is 

s M V ee E ee diag(ei, . . . , e„J , (17) 

since after the transformation from P particles to P 
particles the Hamiltonian in Eq. (fTS"]) has to be diago- 
nal. Multiplying Eq. (jTTJ) from the left by VS' and using 
Eq. (|16a[) yields the eigenvalue equation 

MV = VD , 



where £) = diag(di, . . . ,d ns ) = S'E contains the eigen- 
values of the non-Hermitian matrix M. From Eq. (fTS]) 
below, where we express the Hamiltonian in terms of the 
normal-mode pseudoparticles, it can be seen that the di- 
agonal elements e, correspond to the excitation energies 
of the physical system. Since the energy of the physical 
system must be bounded from below, all have to be 
positive and real, leading to 

e, = dis'i > Vi . 

It will turn out that this stability condition is the only 
point, where the variables s[ of the auxiliary operators 

p** show up. In App.^we show that, if M is completely 
diagonalizablc with real eigenvalues and linear indepen- 
dent eigenvectors, which is of course not generally guar- 
anteed for a non-Hermitian matrix M but necessary from 
the physical viewpoint, then V can be constructed so 
that both requirements of Eqs. (fTBj) and (JTTJ) are fulfilled, 
and we can proceed with our analysis. If M is not com- 
pletely diagonalizable or does not have real eigenvalues, 
the system is unstable, and it favors a different phase, 
which cannot be addressed by the reference system in 
this form. This instability toward a different phase, such 
as superfluidity, has to be cured by extending the refer- 
ence system by proper additional variational parameters, 
as discussed in Sec. IIIII 

In terms of the P operators we obtain for the Hamil- 
tonian 

Pf q = P q 5M q P q = Pi Vl 5M q 7 q P q = PlS'D^ P q 

= ^2 e v(pt,qPv.q ®(Sl,v) +Pu,qPt,q (-^)) 

V 

= E e -^,q^,q- tr 5( £ 'q)- ( 18 ) 

v 

In the last line we have exploited the fact that in order 
for the system to be stable, i. e., the energy be bounded 
from below, all e„ must be positive. 

Inserting this expression in Eq. (TH| yields the Hamilto- 
nian in terms of diagonal normal modes. From this result 
one obtains immediately the grand-canonical ground- 
state energy per cluster 

n = n' + tr 5 (A)--l]rtr ff (P q ) . 

c q 

As discussed, A and P q are diagonal matrices containing 
the poles of the reference Green's function and physical 
Green's function, respectively Therefore, this expression 
being e quiv alent to Eq. (11) in Ref. [ll| (see also Refs. |^, 
[l2l . and l26l) is equivalent to the zero-temperature VCA 
grand potential. 

By using the expression for the Green's function of the 
noninteracting normal modes 

t $a B 

<p«;p« >= , 

M oj — e a 
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we readily obtained the Green's function for the physical 
system 

G q (w) =< a q ; a q >= Q < B q ; B q > Q f 
= QVq « P q ; P q » 7 q +Qt 
= QV^S'uj - S'D^ )- 1 V q t Q t 
= QV^uj - £> q )- 1 1 S'V q t Q t . 

By simple algebra this expression can be rewritten such 
that it is independent of the auxiliary quantities S' and 
V and thus equivalent to Eq. (12) in Ref. [ll] 

GqH = qv^ - y q - 1 Af q Vq)- 1 y- 1 5Qt 

= Q(u - M q )~ W 7 (19) 

where we have used Eqs. (|16[) and (jTTJ). 

We, therefore, succeeded in proving that, for nor- 
mal bosons, the pscudoparticlc approach yields the 
same Green's function and grand potential as VCA. For 
fermions, this was shown in Ref.ll4l see appendix therein. 
This result holds for T = 0, although extension to T > 
is straightforward. 

III. SUPERFLUID PHASE 

When trying to apply VCA to bosonic lattice systems 
in regions of the phase diagram outside the Mott phase, 
one encounters instabilities which manifest in the form 
of noncausal Green's functions, i. e., in spectral functions 
with negative (positive) spectral weight for positive (neg- 
ative) frequencies (J, or in complex poles. Within the 
pscudoparticlc approach these instabilities show up as 
complex eigenvalues or negative diagonal elements of the 
matrix E. This kind of instability is well known in ap- 
proaches based on the bosonic Bogoliubov approxima- 
tion, such as the spin-wave approximation. 

Quite generally, such an instability signals the occur- 
rence of a phase transition toward a new phase. Quite 
often, as in the case of the BH model studied in Sec. IIV1 
the new phase is the superfluid phase, which is accom- 
panied by a Bose-Einstcin condensation. Bose-Einstein 
condensation is described by a finite value of the order 
parameter (or). This suggests to include a source-and- 
drain term in the reference system, which breaks the U (1) 
symmetry of the reference system, leading to the "per- 
turbation" f [sec Eq. 

t = Y, 4 *"(R - R') a R , + $> R f R + 4a R ) , (20) 

where f R = (/i, fa ... /l) t is a vector of size L and is 
identical for all clusters. The index R, however, will be 
kept for notational reasons. 

Due to these terms, the reference system Hamiltonian 
does not conserve particle number anymore. Its eigen- 
states will thus consist of superpositions of states with 



different particle numbers. Numerically, a cutoff in the 
maximum number of boson is necessary in order to solve 
the reference system on the cluster level exactly. We 
again introduce pseudopartiele operators b R connecting 
the ground state with excited states. Note that we can- 
not distinguish between particle or hole states anymore. 
The pseudoparticles are defined by Eq. ([7]) and are con- 
nected to the original boson operators a R by means of 
Eq. |9|). Now, all matrix elements in Eq. (jTUJ) are nonzero 
in general. Therefore, the two sums over [i in Eq. ^ are 
extended to fj. = 1 , . . . , n s , where n s is the number of 
excited states considered in each cluster. 

For the following considerations it is convenient to ex- 
press the boson operators within a Nambu notation. For 
the particle operators we introduce in real space 

which after a Fourier transformation in the cluster vec- 
tors, see Eq. (fT5|) . becomes 

For pseudopartiele operators we have in real space 
Br = (b 1R , b 2R , b na R , b f hR , bl sR f 
and in q space 

Bq — (&l iq ; &2,qi • ■ • ) ^n s ,qi ^1,-q' 1 • • ' ^n s ,-q) ■ 

Similarly to Sec. HH we have an approximate linear 
relation between the A operators and the B operators of 
the form 

a r = qbr + r . 

After the Fourier transformation in the cluster vectors it 
reads 

^ q = Qs q + r q . (21) 

Here, 

r q = v^qr , 

with 

r = (71, 72 ■ ■ • 1L, 7i, 1*2 ■■■ ll ) T , 
and the (2L) x (2n s ) matrix 

Q=( R Z ) 

v \Z* R* ) ■ 

The constants 7, = (V'ol o>% |"0o) j will be nonzero as the 
reference system does not conserve the particle number. 
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In terms of pscudoparticlc operators we can again write 
the reference Hamiltonian for a cluster R, similarly to 
Eq. ((HJ as 



1 



trg(A) 



(22) 



Here, the matrices S and A have a slightly different def- 
inition 



-1) 



n„+l,...,2n s 



S = diag(l,. 

l. 

and 

A = S diag(A£;;, AE' 2 . 



To express the "perturbation" T of Eq. (|20| , we need 
to introduce a similar Nambu notation for the source- 
and-drain terms, which, being R independent, become 
in q space 



AE' AE[, AE> 2 . . . AE' ) 



This gives 

h= N c n' 

N, 



-ytrg(A) --g-trfc- 
/AT 



Av. 



r^r 



2 (F+r + /i.e.) + ^{F^QBo + h.c.) 



±$>t 5A/qBq . 



(26) 



The term linear in B can be eliminated by a shift 



B n = B n + X a 



term of A" q is nonzero. 



where clearly only the q 

Considering only the q = part of Eq. (f2"6")h which we 
term Yq, and plugging in the shifted operators, we obtain 



Y = \{B - X )tSM {B - X ) 



(Nr. 



-{F^Q(B -X ) + h.c.) 



F = 



'N c 5qF 



f 

ftT 



(23) 



After the Fourier transformation in the cluster vectors, 
we can rewrite 



f = f + A = J2(U{f n A l 



1 _ 

2 ^ ^ q 



> Q L q q ' q qJ 

where T q = diag(t q , F q ). 

Replacing the A operators in terms of the B operators 
with the help of Eq. (|2"Tj) . and combining Eq. (|2"2"j) with 

the expression above for T, we finally obtain the com- 
plete Hamiltonian, defined in Eq. ([5]), in terms of pseu- 
doparticles 



H = N c n'+^trg(A)+J2{-l 



tr t r 



ir q T q r q + ^Bl [SA + Q%Q] B. 



SM a 



+ \ [(r q r q + fDqb^ + F q r q + h.c] } . 

The expression can be further simplified by using the fact 
that F and T are equal in all clusters, and thus have only 
q = components. In addition we take advantage of 



Vtrtq = AT c trt(R-R' = 0) = N c trh 



(24) 



The linear term is eliminated by setting 

x = <JW c m^s q^f , 

yielding for the q = term above 



(27) 



1 



Y = -B^SMoBo 



^F f G (0) F. 



where 



G (0) = G q=0 (^ = 0) = -Q M "W 



(28) 



which is the Green's function defined in (|33|) but evalu- 
ated for q = and w = 0. In total we have 

H = N C C+ BlSMqBq (29) 

qGBZ/2 

with the constant terms 

C = n' + \ tr 5 (A) - \ tr h + ±(F*T + h.c.) 

+ ^f T + ^G ia) F . 

In the last term of Eq. (f2"9")l . we restrict the summa- 
tion over half of the Brillouin zone, which we denote by 
q £ BZ/2, and thus removed the factor 1/2 in front of 
the sum. Due to Nambu representation, two summands 
with +q and — q are identical and therefore the restric- 
tion to half of the Brillouin zone is convenient. In our 
convention, the q = term is included in the sum and 
retains the factor 1/2. 



since t(R — R' = 0) = is a pure intercluster term. For 
notational convenience we introduce 



£t = i?t + r t T[) . 



(25) 



A. Condensate density 

Before turning to the diagonalization of the Hamilto- 
nian in Eq. (|2U)) , let us evaluate the condensate density. 
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Since there are no terms linear in B, its expectation value 
(B) vanishes. Therefore, we can immediately calculate 
the condensate density 

(A,} = (A) 

= q + r q = -Qx q + r q 

= y/KS*F + G m (F + T r)], (30) 

where we used Egs. P5]l. (|2"5j) . (|2"7|) and ^25}. We now 
exploit the fact that 

r = (A)' 

is the condensate density in the reference system. From 
the Dyson equation for the Green's function of the phys- 
ical and the reference system we have^ 



G^co)- 1 = G'(u)- 



By multiplying (|30[) with G^ we obtain 



G^ (A) = C;-, 1 (A)' To (A)' + F + f (A)' 



g{-^{aY + f, 



which corresponds to Eq. 



B. Diagonalization of the Hamiltonian 

The Hamiltonian of Eq. (|29p is finally quadratic and its 
diagonalization proceeds in the same way as in Sec. [TXJ 
Again we introduce P operators 

Bq = V q P q , 

and find the solution of the non-Hcrmitian eigenvalue 
equation 

M q v q = y q z? q , 

where Vq satisfies the relation 

VqS'V^S = I . 

The diagonal matrix S', which is in principle q-dependent 
as well, consists of +1 or — 1 terms. It is chosen accord- 
ing to the prescription derived in App. [3] The stability 
condition is again that the pseudoparticle eigenenergics 

S'Dq = diag(e iq , . . . , e 2 „ sq ) 

are all positive. The physical Hamiltonian in terms of 
P-particles now reads 

k = E PlS'D^ + NcC 

qGBZ/2 

= E Ew4^,q- E 9(D^)+N C C. 
qGBZ/2 V qGBZ/2 

(31) 



From that wc readily obtain (see App. [Bj the grand po- 
tential per cluster of the physical system f2, which is the 
ground state expectation value (H) /N c 



tr<?(A) 



1 

~Nr. 



E 

qGBZ/2 



trh 



+ i(At)G^(A)-i(At)'G' ( - ) 1 <A)' 



(32) 



By considering the fact that A and Z) q contain the 
poles of G' and G, respectively, we conclude that, in the 
T -> limits 



lim 

T-s-0 



lTrln(-G')-^Trln(-G) 



A, 



tr g(A) 



E 9(D«) 

qGBZ/2 



Thus, Eq. (|32j) is equivalent to Eq. ([T]) in the introduction 
in the T = limit. An extension to T > is straight 
forward. 

The connected Green's function now contains anoma- 
lous contributions, but formally is obtained as in Eq. (fT9|) . 

G q M =« A q ; A\ » c = Q « P q ; P q » 

= QV q « P q ; Pi » y q tQt 

= QK 1 (5' W -5'P q )- 1 y q t Q t 

= 0K 1 ( w -p q )- 1 y q - 1 5Q t 

= Q(u - M q )- W , (33) 

where we have neglected the shifts T and X since 
they only contribute to disconnected parts. Notice that 
Eq. p3|) is a 2L x 2L matrix in Nambu and cluster-site 
space. The q vectors above refer to the reduced Bril- 
louin zone originating from the cluster tiling, therefore 
G is expressed in a mixed representation. In translation- 
invariant systems, the Green's function is expected to be 
diagonal in the wave vectors k of the full Brillouin zone. 
This symmetry is notoriously broken in cluster methods 
such as VCA or C-DMFT. In order to obtain a k-diagonal 
2x2 Nambu Green's function G(k, u>) we need to apply 
a periodization prescription^ This gives 



G(k, w) = vtG k (w)v 



k J 



where 



-i k ri 



o 



-ikr_L 







-tkri 



o 

-ikri 



(34) 

and Vi is the position of site i within the cluster. 

A nontrivial test for VCA is the noninteracting limit, 
for which this approximation becomes exact. In Ap- 
pendix [C] we carry out this check for the noninteracting 
BH model, i.e., we set U = 0, and for a reference sys- 
tem consisting of single-site clusters. In this test case the 
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grand potential f2 of the physical system can be evaluated 
analytically both using the VCA prescription as well as 
directly from the Hamiltonian of noninteracting lattice 
bosons. 



C. Particle density and momentum distribution 

The total particle density is defined as 

q i 

where N = N c L is the total number of lattice sites 
present in the physical system. The particle density can 
be easily expressed in Nambu formalism 

11= JN EE^qO + K-q4-q)) _ \ 
q i 

= ~2 + 2N^ A ^ 
q 

= ~\ + ^ E^qWQK^q) + (4,)* <4,» 

q 

= ~\ + ^ E tr[e(-i? q )V q tQt gy q ] + J_ (^) (A) , 

(35) 

where the last term describes the contribution from the 
condensate, which can be deduced from Eq. (f5D|) . The 
term with the sum over q can be rewritten to obtain the 
known form of the particle density^ 

n = -\ 2^E tr t (- i5 q) 5 % t Q t ^q] + 51 ^ ^ 
q 

= ~\- E tr [ (- i? q)^ 1 ' 5£ 3 t£ 3 V q] + ^ ^ ^) 

q 

The momentum distribution n(k) can be extracted by 
the Fourier transform within the cluster leading to 

n(k) = — — + — (A*) (A) 

+ ^tr[v] £ gF k e(-^ k )v k t g t v k ], 

where v£ is given by Eq. ([34]) . 

IV. APPLICATION TO THE BOSE-HUBBARD 
MODEL 

In this section, we present the first nontrivial applica- 
tion of the extended VCA theory to the two-dimensional 
BH (BH) model and compare the results with unbiased 
quantum Monte Carlo (QMC) calculations. The BH 




0.02 0.04 0.06 0.08 

t/U 



FIG. 1. (Color online) Phase boundary for the first three 
Mott lobes corresponding to filling n = 1, 2 and 3. The data 
for the first two Mott lobes have been published in Ref. [Til . 
Static quantities are evaluated along the dashed line, i.e., 
for t/U — 0.02 and fi/U ranging from to 3, whereas, the 
dynamic single-particle spectral function is evaluated at t = 
0.07 and /j, — 0.4, see mark x. 

Hamiltonian , 15 ! 18 which describes strongly correlated lat- 
tice bosons, reads 

^ = -* E a » a i + ^E hi - !) - mE ft * ' 

where a\ (a i ) creates (destroys) a bosonic particle and 
hi = a\ a i counts the number of particles at lattice site i. 
The parameter t is the hopping strength, which originates 
from the overlap of the localized wave functions belonging 
to lattice sites i and j, respectively. The first sum (indi- 
cated by angle brackets) is restricted to ordered pairs of 
nearest neighbor sites. The repulsive on-site interaction 
is termed U, and \x is the chemical potential, which con- 
trols the particle number. For increasing ratio t/U the 
system undergoes a quantum phase transition from the 
Mott to the supcrfluid phase. We evaluate static quan- 
tities, such as the particle density n and the condensate 
density n c as well as the dynamic single-particle spec- 
tral function A(k, ui). The phase boundary of the first 
three Mott lobes as obtained in VCA is shown in Fig.[T] 
The data for the first two lobes have been published in 
Ref. [TJJ. Static quantities are evaluated for constant hop- 
ping strength t/U = 0.02 and distinct values of the chem- 
ical potential fj,/U ranging from to 3, scanning through 
various Mott lobes separated by the supcrfluid phase; 
see the dashed line in Fig.[T] The single-particle spec- 
tral function is evaluated for the parameter set marked 
by x in Fig.[TJ which is located in the supcrfluid phase 
close to the tip of the first Mott lobe. For the numeri- 
cal evaluation we used the chemical potential fi' and the 
strength of the source-and-drain coupling term F of the 
reference system as variational parameters. If not stated 
differently, the reference system consists of a cluster of 
size L = 2 x 2. 
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FIG. 2. (Color online) Total particle density n, condensate 
density n c , and density of the particles which are not con- 
densed n — n c evaluated along the dashed line shown in Fig.[T] 
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FIG. 4. (Color online) Single-particle spectral function 
A(k, w) evaluated at t/Z7 = 0.07 and /Lt/f7 = 0.4. The colored 
density plot corresponds to VCA results and the dots with 
errorbars to latest QMC results of Ref. I32I 
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FIG. 3. (Color online) Comparison of the total particle den- 
sity n evaluated by means of VCA and QMC for parame- 
ters along the dashed line in Fig.[T] (t/U = 0.02). The inset 
compares VCA and QMC results for the condensate fraction 
n c /n. VCA results are obtained for reference systems of size 
L = 1 X 1 and L = 2 x 2 and essentially infinitely large phys- 
ical systems. QMC results are obtained for physical systems 
of size 32 x 32 inverse temperature U/T = 128. 



The total particle density n evaluated using Eq. ([55)1 
is shown in Fig. [5] along with the condensate density 
n c = (A') (A) /2L, and the density of the particles which 
are not condensed n — n c . From Fig.[3]it can be observed 
that the particle density n evaluated for reference systems 
of size L = 1 X 1 and of size L = 2 x 2 are almost iden- 
tical. The same holds for the condensate fraction n c /n, 
which is shown in the inset of Fig-El In the same figure, 
we also compare our results with QMC calculations. The 
densities obtained from the two methods show an excel- 
lent agreement. The QMC data have been obtained for a 
system of size 32 x 32 and temperature U/T = 128 using 
the ALPS library^ , and the ALPS applications .21 

The single-particle spectral function A(k, u>) evaluated 



for the parameter set. marked by x in Fig.[IJ i. e., in the 
supcrfluid phase close to the tip of the first Mott lobe, 
is depicted in Fig. [4] The colored density plot corre- 
sponds to VCA results and the dots with errorbars to 
latest QMC results of Ref. HI The VCA spectral func- 
tion A(k, ui) consists of four bands, which is in agreement 
with results obtained by means of a variational mean field 
calculation^ a strong coupling approach^ and random 
phase approximation (RPA) calculations . 35 ! 36 The ad- 
vantage of VCA in comparison to the above mentioned 
approaches is that the results can be systematically im- 
proved by increasing the cluster size of the reference sys- 
tem. For each wave vector k the weight is concentrated 
in one of the two bands present at positive and nega- 
tive energy, respectively. We observe that the outer two 
modes exhibit a wide gap at k = 0, which is approxi- 
mately of size U. The inner two, low-energy modes are 
also gapped at k = 0. However, the gap is tiny, and 
away from k = the spectrum quickly develops a linear 
behavior, which is in agreement with the expected dis- 
persion of Goldstonc modes. The failure in obtaining a 
gapless long-wavelength excitation is a common problem 
of conserving approximations, i. e., of approximations for 
which macroscopic conservation laws are fulfilled. Simi- 
lar aspects occur in dynamical mean-field theory calcu- 
lations of two-component ultracold atoms as welli^i In 
VCA there exists the additional possibility to systemati- 
cally improve the obtained results by increasing the clus- 
ter size L of the reference system. Figure [S] compares the 
k = gap of the inner modes for reference systems of size 
L = 1 X 1 and L = 2 x 2. The gap is evaluated along the 
dashed line shown in Fig.[TJ The first observation is that 
the gap present in the condensed phase is almost an order 
of magnitude smaller than the gap in the Mott phase. It 
vanishes at the Mott-to-supcrfluid transition and, most 
importantly, shrinks with increasing cluster size L. This 
behavior signals convergence toward the correct result. 
In Fig.2] we also compare our VCA results for the 
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FIG. 5. (Color online) Gap of the inner modes present in 
the single-particle spectral function measured at k = and 
evaluated along the dashed line shown in Fig.[TJfor reference 
systems of size L = 1 x 1 and L — 2 x 2, respectively. 



single-particle spectral function to latest QMC results 
obtained in Ref. [H. In this figure QMC results are in- 
dicated by dots with errorbars, which quantify the peak 
position of the spectral weight. Overall, we find good 
agreement in the low-energy spectrum. Only very close 
to k = the two results differ slightly and the QMC dis- 
persion possesses the correct gapless behavior. The QMC 
spectral function, exhibits only two instead of four bands. 
This is, however, not surprising since for the considered 
parameter set and at a specific wave vector k the weight 
of one positive (negative) energy band dominates dras- 
tically over the other one located at positive (negative) 
energy. Thus the four bands are extremely difficult to re- 
solve by means of the maximum entropy method, which 
has been used to infer the spectra from QMC data; see 
Ref. HH for details concerning the QMC results. This 
reference also contains a comparison between VCA data 
and QMC data for the spectral function evaluated in the 
Mott phase, where the results obtained from the two ap- 
proaches coincide very well for all k values. 

We also evaluated the particle density n for the pa- 
rameter set in the superfluid phase used in Fig. 2] and 
compared it to the QMC results. VCA yields n = 1.0321 
in excellent agreement with the QMC result n^ MC = 
1.03068(2) obtained at U/T = 128 for a system of size 
32 x 32. 

In the following we provide some additional remarks 
on the fact that the single-particle spectral function, 
obtained within our approach, is gapped in the long 
wavelength limit, i.e., close to k = 0, for modes 
which ought to be identified with the Goldstone modes. 
This issue is a commonly known problem of conserving 
approximations^ One condition for an approximation 
to be conserving is, for example, to be ^-derivable and 
self-consistent (see Ref. I39rl43l for details). VCA is $- 
derivable but not self-consistent: The self-energy is the 
derivative of a functional of the Green's function, but the 



latter is not the Green's function obtained from Dyson's 
equation. Thus VCA is not completely conserving. How- 
ever, many conservation laws are fulfilled at the sta- 
tionary point of the self-energy functional, depending on 
which variational parameters are taken into account (see 
Ref. for a more detailed discussion). 



To obtain a gapless spectrum, a system of condensed 
bosons has to fulfill an independent condition, which 
is the Hugenholtz-Pincs theorem i 42 i 45 i 46 There are only 
very few systematic approximation schemes which sat- 
isfy both conditions simultaneously. One notable excep- 
tion occurs for interacting bosons composed of paired 
fermions. In this case, a consistent and gapless ap- 
proximation can be developed provided the theory is ex- 
pressed in terms of the constituent fermions 4£ In a dif- 
ferent work 4 ^ it was suggested to include an additional 
Lagrange multiplier in the form of a chemical potential, 
in order to explicitly enforce the Hugenholtz-Pines con- 
dition. Unfortunately, the Hugenholtz-Pines theorem is 
not fulfilled in VCA, and thus the low-energy modes of 
the single- particle spectral function are gapped in the 
long wavelength limit. Yet, the gap present in the VCA 
single-particle spectral function is small, and the spec- 
trum quickly develops a linear behavior reminiscent of 
the gapless and linear Goldstone modes. Furthermore, 
in VCA there exists the possibility to systematically im- 
prove the results by increasing the cluster size of the ref- 
erence system. 



It is also interesting to mention that the related strong 
coupling approximation RPA, which yields a gapless sec- 



trum, yet is not conservin 



can be obtained within 



certain limits of the extended VCA formalism. Specif- 
ically, the limits to consider arc (i) to use clusters of 
size L = lxl, (ii) not to use the chemical potential 
/i as variational parameter and (iii) to determine the 
source-and-drain coupling strength F self-consistently 
within a mean-field approach, whereby intercluster hop- 
ping terms a\ ctj are replaced with their mean-field value 



at) a,j+a 



in the reference Hamiltonian. This leads 



to the selfconsistency condition F = zt (A) , where (A) 
is given by Eq. (|30p and z is the coordination number of 
the lattice. Our formalism provides a natural way to im- 
prove on RPA in a gapless, yet nonconserving, way by 
simply increasing the cluster size L and fixing F using 
the mean-field condition discussed above. However, it 
has to be emphasized that VCA yields much better re- 
sults than RPA, even if RPA is extended to clusters of 
size L. Specifically, the particle density, the condensate 
density and the location of the phase boundar y 1 ^ 49 can 
be determined much more accurately by means of VCA 
only because we allowed for a variation in the chemical 
potential fx', i.e., allowed for macroscopic conservation 
laws to be fullfilled. 
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V. CONCLUSIONS 

In the present paper, we introduced a pseudoparticle 
formalism for interacting bosonic systems, and showed 
that the results of the variational cluster approach can 
be derived within this formalism. We used it to extend 
the variational cluster approach to the superfiuid phase 
of strongly correlated lattice bosons. We derived expres- 
sions for the grand potential and for other quantities, 
which are necessary to investigate the superfiuid proper- 
ties. Our results suggest that the pseudoparticle formal- 
ism is a quite versatile approach, which can be applied 
to a large variety of other problems. 

As a first nontrivial application of the extended version 
of the variational cluster approach we choose the two- 
dimensional Bose-Hubbard model and evaluated static 
quantities such as the total particle density and the con- 
densate density, as well as the dynamic single-particle 
spectral function. We compared the single-particle spec- 
tral function with recent Quantum Monte-Carlo results 3 -^ 
and found good agreement between the two approaches. 
It has to be pointed out that our extended variational 
cluster approach, while fulfilling many conservation laws, 
does not fulfill the Hugenholtz-Pines theorem. From this 
fact follows that the low-energy excitations of the spec- 
trum have a small but nonzero gap in the long wavelength 
limit. This is a common aspect, which is already present 
in theories of the dilute Bose gas i 38 ' 50 ' 51 However, for 
wavevectors away from k = the spectra obtained within 
this approach quite soon exhibit a correct linear behavior 
and agree very well with the Quantum Monte- Carlo re- 
sults. Moreover, the gap shrinks with increasing cluster 
size, corroborating that the variational cluster approach 
becomes exact in the infinite cluster limit. Due to the 
fact this approach fulfills several conservation laws, the 
particle density, the condensate density as well as the 
phase boundar y 11 ' 49 delimiting the Mott from the su- 
perfiuid phase can be evaluated very accurately. In the 
present paper we demonstrated, that our variational clus- 
ter approach results for the densities evaluated in both, 
the Mott and the superfiuid phase, match perfectly with 
Quantum Monte-Carlo results. 
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Appendix A: Procedure to construct V and S' 

Here, we outline how the two conditions on V given 
in Eqs. (|16[) and (JTTJ) can be achieved and how S' can be 



constructed. We start out from the eigenvalue equation 
for the non-Hcrmitian matrix M 

AIV = VD . 

As already argued in Sec. UH from the physical view- 
point we can only proceed if the eigenvector-matrix V is 
nonsingular and if all eigenvalues are real, as the system 
would otherwise be unstable. Hence we can express the 
Hcrmitian diagonal matrix of eigenvalues as 

D = V~ X MV . 

The first condition of Eq. (|16[) requires that the Hermitian 
matrix 

X = V^SV 

be diagonal with diagonal elements Xu = ±1. Multiply- 
ing the two Hermitian matrices and exploiting the Her- 
miticity of SM results in 

XD = V^SMV = (XDy = DX ; =S> [X, D] = . 

Commuting Hcrmitian matrices have a common set of 
orthonormal eigenvectors. The matrix D is already di- 
agonal. Hence for indices belonging to nondegenerate 
eigenvalues, X is also diagonal. Within the set of indices 
belonging to a degenerate eigenvalue, the corresponding 
Hcrmitian submatrix of X can be diagonalized by a uni- 
tary transformation U . In the following we term the di- 
agonalized matrix as X' . The diagonalization also results 
in a new matrix V = VU of eigenvectors. We still have 
V~ 1 MV = D, but now 

V^SV = X' = diag(x' 1 ,...,x' L ) (Al) 

V f SMV = E' = X'D = diag(a;i<2i, . . . , x' ns d ns ) . 

For the condition Eq. (fTB|) we still need to ensure that 
x' a = ±1. Provided no x' a vanishes, which we will show 
below, this can easily be achieved by a suitable normal- 
ization of the column vector of V — > V = VZ , wi th Z 
being a diagonal matrix, defined as Z aa = 1/ We 
eventually have 

V~ X MV = D = diag(di ,...,d L ) 

V^SV = Z^X'Z = S' = diag(s' 1 , . ..,s' L ) 
V^SMV = E = diag(ei, . . . ,e„J . 

We are merely left with the proof that 

x' a = v a, Sv" ^ , (A2) 

where v Q stands for the ath column of V . To this end 
we assume ad absurdum that v°'Sv° =0. In this case, 
v Q would belong to the (n s — l)-dimensional space S a 
orthogonal to the vector Sv a . According to Eq. (|A1[) the 
vectors v , ... v"" 1 , v Q+1 , . . . v" s also belong to S a and 
they are linear independent. Thus they span S a . Due to 
the fact that all vectors v 1 , ... v' ls are linear indepen- 
dent, v a cannot belong to S a , which proves Eq. (|A2|) . 
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Appendix B: Grand potential 



Comparison reveals 



In this appendix we derive Eq. (|32|) . Starting out from 
Eq. ([21]) we get 



C qeBZ/2 

= n' + itr fl (A)-^- £ flW + ^F+r + ^c.) 

(Bl) 



qGBZ/2 



_ i tr/l+ i r tfor+iFtG ( o)F. 

We now evaluate the quantity 

W = (A^)G^(A)-(A^G' ( -f(Ay 

= r t G (0 5r + r t i= :, + F t r+ 

= (T\F + f T) + h.c.) + F^G {0) F - T^foT 
= (T^F + h.c.) + T^f T + FG (0) F . 

Comparison with (|B1[) gives 
= fi' + itr 5 (A)-i- 9{D«) + ±W-l 



trh 



qSBZ/2 



which is the expression for the grand potential stated in 
Eq.®. 



a = —fj, 

x = -f/a = f/(jf 
c=-a\x\ 2 = \f\ 2 /p>. 

The Hamiltonian H', rewritten by means of the shifted 
operators, is given by 

H' = -n'tfa + \f\ a /»'. 

As discussed before we choose // < 0. The eigenencrgics 
obtained form the Schrodinger equation are 

it'\P) = (-»'i>+\f\ 2 /p')\i>)=E' u \i>) . 

For negative chemical potential p! the ground state is 
|V'o) = 1 6) and its energy Eq = \f\ 2 /p'. The eigenstates 
of H' are number states, therefore the shifted creation 
and annihilation operators act on them in the usual way 

a \v) = \v — 1) 



\r, 



\S) = VP + T\v + 1) . 

To evaluate the Q matrices we apply the original op- 
erators a on the eigenstates of H' 

a \v) = (5 - f/p') \v) =yfi\v-l)- f/p' \v) . 
With that we obtain 



Appendix C: Zero-interaction limit 

The zero-interaction limit turns out to be a nontrivial 
check for VCA. For U = 0, the BH model can be solved 
analytically as it reduces to 



H = -t ^ a\a j -pJ2- 



The chemical potential p has to be smaller than — It in 
order to prevent infinitely many particles in the ground 
state. Taking this into account, the grand potential at 
zero temperature is O = 0. In the zero- interaction limit 
VCA/CPT yields exact results. Thus the pseudoparti- 
cle formalism can be checked by applying this limit. For 
reference systems W which consist of a single site the 
calculations can be done analytically. Under these con- 
siderations the Hamiltonian H' reads 

H' = -p! o f a - {a) f + f* a) . 

It can be solved by introducing shifted operators a = a+x 
and by "completing the square" 

it' = -p' a) a - (a f / + /*a) 

= aa) a + c — a{a) + x*)(a + x) + c 

= aa' a + a{a)x + x*a) + a \x\ 2 + c . 



<0| o |£>> = 1 
(z>|a|6) = 
(0| a |0) = -////■ 
Writing down the expressions in matrix form yields 



A = s( E >- E ° „ 



E[ - E 



-p' 

p! 



Using the expressions above and the relation A = QB + T 
we obtain for the pseudoparticle operators 

B = Q-\A-r) = A . 

Next, we evaluate the grand potential from Eq. (|B1[) . 
where we obtain 

fi = fi' + i t r 5 (A)-i- ]T g(Dj + l(Fir + h.c) 



qGBZ/2 



- tr h+ ^f T - -F f Q M^SQpF 
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by employing Eq. (f28|) . We calculate parts A-F of ft sep- 
arately 



A: fi' + -tr< 7 (A) = |/| 2 / /t / + ///2 



qGBZ/2 



C: -(F^ + h.c.) = -2\f\ 2 /»> 

D: itr/i=( M / - M )/2 

E: irtr r = |/| V - ^ - 2;)/ M ' 2 

F: i^QM^ = (At + 2 i)/ / / 2 . 



where we used that 



t 



and i q = iT q = // — /U — 2t cos q. Since M q is already 
diagonal we can readily evaluate part B as sum over the 
negative eigenvalues, which is (i + 2icosq, since y, < 
—2t. When summing over half of the q values the second 
term of the eigenvalue containing cos q is zero. For the 
calculation of part F we need the inverse of Mo , which is 
simply 



M„ 



and F. which reads 







fj,+2t 



In order to evaluate part B we need the matrix M q , 
which is given by 



Ma = A + SQ^f^Q 



—fj, — 2icosq 

fi + 2t cos q J ' 



F = F + T T = + 2t)/n' 



f* 



Collecting all terms yields the grand potential O = 0, 
which is identical to the result obtained from the direct 
calculation. 
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